package scu.maqiang.fes;

import java.util.Arrays;
import java.util.HashSet;
import java.util.function.Function;

import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarCFunc;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.*;
import scu.maqiang.peridynamics.interpolation.PDApproximation;

public class FES1L21 extends FES  {

	public FES1L21(Mesh1L2 mesh1D) {
		super(mesh1D, 1);
		elementElasticStrain = new double[mesh1D.getNt()];
		elementPlasticStrain = new double[mesh1D.getNt()];
		elementStress = new double[mesh1D.getNt()];
		eleStiff = FES1L21::elementStiff;
		eleMass = FES1L21::elementMass;
		eleSource = FES1L21::elementSource;
		eleFlux = FES1L21::elementFlux;
	}
	
	/**
	 * ����д��һ����ͨ�÷���
	 * @param u
	 * @param x
	 * @param eleIdx
	 * @param eleVal
	 * @return
	 */
	public double getValue(double[] u, double[] x, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
	    double[] uu = new double[2];
	    double value = 0;
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coordRegion);	    	
	    	if (CG.pointIn1DL2(coordRegion, x, eleVal)) {
	    		for(int j = 0; j < 2; j++) {
	    			uu[j] = u[ele[j]];
	    			eleIdx[j]= ele[j];
	    			value += uu[j] * eleVal[j];
	    		}
	    		return value;
	    	}
	    }
	    throw new IllegalArgumentException("The Point ( " + Arrays.toString(x) + " ) is not in the domain!");
	}
	
	
	public void getValue(double[] u, double[] xyz, DoubleIntegerPair data, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		double[][] coord = new double[2][1];
		int[] ele = mesh.getElement(data.second);
		mesh.getCoordPerElement(ele, coord);
		if (CG.pointIn1DL2(coord, xyz, eleVal)) {
			for (int j = 0; j < 2; j++) {
				// eleIdx[j] = ele.idx[j];
				value += u[ele[j]] * eleVal[j];
			}
			data.first = value;
			return;
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coord);
				if (CG.pointIn1DL2(coord, xyz, eleVal)) {
					for (int j = 0; j < 2; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					data.first = value;
					data.second = i;
					return;
				}
			}

		}
		throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}
	
//	public void assembleStiff(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A) {
//		assembleGlobalMatrix(coefFunc, param, FES1L21::elementStiff, tp, A);
//	}
//	public void assembleStiff(double[] constCoef, BVPType tp, SRMatrix A) {
//		assembleGlobalMatrix(constCoef, FES1L21::elementStiff, tp, A);
//	}
//
//	public void assembleStiff(double[][] vecCoef, BVPType tp, SRMatrix A) {
//		assembleGlobalMatrix(vecCoef, FES1L21::elementStiff, tp, A);
//	}

//	public void assembleStiff(double[] constCoef, double[][] vecCoef, BVPType tp, SRMatrix A) {
//		assembleGlobalMatrix(constCoef, vecCoef, FES1L21::elementStiff, tp, A);
//	}

	public static void elementStiff(int[] ele, double[][] xx, double[][] kxx, BVPType tp, double[][] eleMatrix) {
		double[] x = MVO.col(xx, 0);
		double L = Math.abs(x[0] - x[1]);
		double xi, detJ;
		double[] h = {0.0, 0.0};
		double[] dx = {0.0, 0.0};
		double coef, c, wxarea;
		double[] coef0, coefX, coefX2, coefX3, coefU, coefU2, coefU3, coefdU, coefdU2, coefdU3, u;
		MVO.fill(eleMatrix, 0.0);
		switch (tp) {
			case COMMON:
				coef = 0.5 * MVO.sum(kxx[0]) / L;
				eleMatrix[0][0] = eleMatrix[1][1] = coef;
				eleMatrix[0][1] = eleMatrix[1][0] = -coef;
				break;
			case AXISYMMETRIC:
				coef = (MVO.sum(kxx[0]) * MVO.sum(x) + MVO.dot_product(kxx[0], x)) / (6.0 * L);
				eleMatrix[0][0] = eleMatrix[1][1] = coef;
				eleMatrix[0][1] = eleMatrix[1][0] = -coef;
				break;
			case SPHERICAL:
			case NONLINEAR_2:
				//系数形式为 a + b x + c x^2 + d x^3 + e u + f u^2 + g u^3 + h du/dx + k (du/dx)^2 + l (du/dx)^3 的非线性切线刚度矩阵
				coefU   = kxx[4];
				coefU2  = kxx[5];
				coefU3  = kxx[6];
				coefdU  = kxx[7];
				coefdU2 = kxx[8];
				coefdU3 = kxx[9];
				u       = kxx[10];
				for(int i = 0; i < 4; i++) {
					xi = Constants.GP4[i];
					detJ = CG.ISOMap1D(CG::L2ShapeFunction, x, xi, h, dx);
					double c1 = 0.0;
					double uh = MVO.dot_product(u, h);
					double dxUh = MVO.dot_product(dx, u);
					c1 += MVO.dot_product(coefU, h) * dxUh;
					c1 += 2 * MVO.dot_product(coefU2, h) * uh * dxUh;
					c1 += 3 * MVO.dot_product(coefU3, h) * uh * uh * dxUh;

					double c2 = 0.0;
					c2 += MVO.dot_product(coefdU, h) * dxUh;
					c2 += 2 * MVO.dot_product(coefdU2, h) * dxUh * dxUh;
					c2 += 3 * MVO.dot_product(coefdU3, h) * dxUh * dxUh * dxUh;
					wxarea = Constants.GW4[i] * detJ;
					for(int si = 0; si < 2; si++) {
						for(int sj = 0; sj < 2; sj++) {
							eleMatrix[si][sj] += dx[si] * h[sj] * wxarea * c1 + dx[si] * dx[sj] * wxarea * c2;
						}
					}
				}
			case NONLINEAR_1:
				//系数形式为 a + b x + c x^2 + d x^3 + e u + f u^2 + g u^3 + h du/dx + k (du/dx)^2 + l (du/dx)^3 的非线性直接迭代矩阵
				coef0   = kxx[0];
				coefX   = kxx[1];
				coefX2  = kxx[2];
				coefX3  = kxx[3];
				coefU   = kxx[4];
				coefU2  = kxx[5];
				coefU3  = kxx[6];
				coefdU  = kxx[7];
				coefdU2 = kxx[8];
				coefdU3 = kxx[9];
				u       = kxx[10];
//				double avA = MVO.average(coef0) / L;
//				double avBx = (2 * coefX[0] * x[0] + coefX[0] * x[1] + coefX[1] * x[0] + 2 * coefX[1] * x[1]) / (6 * L);
//				double avCx2 = 0.25 * (coefX2[0] * x[0] * x[0] + coefX2[1] * x[1] * x[1]) / L
//						      + (coefX2[0] + coefX2[1]) * x[0] * x[1] / (6 * L)
//						      + (coefX2[0] * x[1] * x[1] + coefX2[1] * x[0] * x[0]) / (12 * L);
//				double avDx3 = ( (4 * coefX3[0] + coefX3[1]) * x[0] * x[0] * x[0]
//				               + (3 * coefX3[0] + 2 * coefX3[1]) * x[0] * x[0] * x[1]
//				               + (2 * coefX3[0] + 3 * coefX3[1]) * x[0] * x[1] * x[1]
//				               + (coefX3[0] + 4 * coefX3[1]) * x[1] * x[1] * x[1]) / (20 * L);
//				double avEu = (2 * coefU[0] * u[0] + coefU[0] * u[1] + coefU[1] * u[0] + 2 * coefU[1] * u[1]) / (6 * L);
//				double avFu2 = 0.25 * (coefU2[0] * u[0] * u[0] + coefU2[1] * u[1] * u[1]) / L
//						      + (coefU2[0] + coefU2[1]) * u[0] * u[1] / (6 * L)
//						      + (coefU2[0] * u[1] * u[1] + coefU2[1] * u[0] * u[0]) / (12 * L);
//				double avGu3 = ((4 * coefU3[0] + coefU3[1]) * u[0] * u[0] * u[0]
//						      + (3 * coefU3[0] + 2 * coefU3[1]) * u[0] * u[0] * u[1]
//						      + (2 * coefU3[0] + 3 * coefU3[1]) * u[0] * u[1] * u[1]
//						      + (coefU3[0] + 4 * coefU3[1]) * u[1] * u[1] * u[1]) / (20 * L);
//				dxU = (u[1] - u[0]) / L;
//				double avHDu  = 0.5 * (coefdU[0] + coefdU[1]) * dxU / L;
//				double avKDu2 = 0.5 * (coefdU2[0] + coefdU2[1]) * dxU * dxU / L;
//				double avLDu3 = 0.5 * (coefdU3[0] + coefdU3[1]) * dxU * dxU * dxU / L;
//				double totalCoef = avA + avBx + avCx2 + avDx3 + avEu + avFu2 + avGu3 + avHDu + avKDu2 + avLDu3;
//				eleMatrix[0][0] += totalCoef;
//				eleMatrix[0][1] += - totalCoef;
//				eleMatrix[1][0] += - totalCoef;
//				eleMatrix[1][1] += totalCoef;
				for(int i = 0; i < 4; i++) {
					xi = Constants.GP4[i];
					detJ = CG.ISOMap1D(CG::L2ShapeFunction, x, xi, h, dx);
					c = MVO.dot_product(coef0, h);
					double xh = MVO.dot_product(x, h);
					c += MVO.dot_product(coefX, h) * xh;
					c += MVO.dot_product(coefX2, h) * xh * xh;
					c += MVO.dot_product(coefX3, h) * xh * xh * xh;
					double uh = MVO.dot_product(u, h);
					c += MVO.dot_product(coefU, h) * uh;
					c += MVO.dot_product(coefU2, h) * uh * uh;
					c += MVO.dot_product(coefU3, h) * uh * uh * uh;
					double dxUh = MVO.dot_product(u, dx);
					c += MVO.dot_product(coefdU, h) * dxUh;
					c += MVO.dot_product(coefdU2, h) * dxUh * dxUh;
					c += MVO.dot_product(coefdU3, h) * dxUh * dxUh * dxUh;
					wxarea = Constants.GW4[i] * detJ * c;
					for(int si = 0; si < 2; si++) {
						for(int sj = 0; sj < 2; sj++) {
							eleMatrix[si][sj] += dx[si] * dx[sj] * wxarea;
						}
					}
				}
//				System.out.println(MVO.toString(eleMatrix));
//				System.exit(1);
			case NONLINEAR_POLYNORMIAL_2:

			case NONLINEAR_POLYNORMIAL_1:
				break;
			default:
				throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
	public void assembleElasticStiff(ScalarFunc[] ENuFunc, double[] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(ENuFunc, param, FES1L21::elementElasticStiff, tp, A);
	}
	
	public void assembleElasticStiff(double[] ENu, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(ENu, FES1L21::elementElasticStiff, tp, A);
	}

	public void assembleElasticStiff(double[][] vecCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(vecCoef, FES1L21::elementElasticStiff, tp, A);
	}

	public static void elementElasticStiff(int[] ele, double[][] xx, double[][] ENu, BVPType tp, double[][] eleMatrix) {
	    double[] E = ENu[0];
	    double[] Nu = ENu[1];
	    double[] x = MVO.col(xx, 0);
	    double len = Math.abs(x[1] - x[0]);
	    double mE = MVO.average(E);
	    double mNu = MVO.average(Nu);
	    MVO.fill(eleMatrix, 0.0);
	    double[][] B = new double[2][2];
		switch (tp) {
			case COMMON -> {
				double lambda = mE * (1 - mNu) / ((1 - 2 * mNu) * (1 + mNu));
				double coef = lambda / len;
				eleMatrix[0][0] = eleMatrix[1][1] = coef;
				eleMatrix[0][1] = eleMatrix[1][0] = -coef;
			}
			case AXIS_PLANESTRESS, AXIS_PLANESTRAIN -> {
				double[][] D = Materials.ElasticConstituteMatrix(mE, mNu, tp);
				double L0 = MVO.sum(x) * 0.50;
				B[0][0] = -1.0 / len;
				B[0][1] = -B[0][0];
				B[1][0] = B[1][1] = 0.5 / L0;
				MVO.mul(MVO.matmul(MVO.transpose(B), MVO.matmul(D, B)), len * L0, eleMatrix);
			}
			case SPHERICAL -> {
				double[][] DD = Materials.ElasticConstituteMatrix(mE, mNu, tp);
				double L00 = MVO.sum(x) * 0.50;
				B[0][0] = -1.0 / len;
				B[0][1] = -B[0][0];
				B[1][0] = B[1][1] = 0.5 / L00;
				double squareX = (x[0] * x[0] + x[1] * x[1] + x[0] * x[1]) / 3;
				MVO.mul(MVO.matmul(MVO.transpose(B), MVO.matmul(DD, B)), len * squareX, eleMatrix);
			}
			default -> throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
	public void assembleHeatThetaScheme(ScalarFunc[] coefFunc, double[] param, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		double dt = param[0];
		double theta = param[1];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int label = mesh.getElementLabel(i);
			for (int j = 0, nj = nPerEle; j < nj; j++) {
				kxx[0][j] = coefFunc[0].action(coordRegion[j], label, param);
				rho[0][j] = coefFunc[1].action(coordRegion[j], label, param);
			}
			elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
			for(int ii = 0; ii < dofPerElement; ii++ ) {
				for(int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}
	
	public void assembleHeatThetaScheme(double[] coefFunc, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		Arrays.fill(kxx[0], coefFunc[0]);
		Arrays.fill(rho[0], coefFunc[1]);
		double dt = coefFunc[2];
		double theta = coefFunc[3];		
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
//			System.out.println(MVO.toString(eleMatrixMass));
//			System.out.println(MVO.toString(eleMatrixStiff));
//			System.exit(0);
			for(int ii = 0; ii < dofPerElement; ii++ ) {
				for(int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}

	public void assembleSource(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS) {
		assembleGlobalVector(coefFunc, param, FES1L21::elementSource, tp, RHS);
	}
	
	public void assembleSource(double[] constCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(constCoef, FES1L21::elementSource, tp, RHS);
	}

	public void assembleSource(double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(vecCoef, FES1L21::elementSource, tp, RHS);
	}

	public void assembleCSource(ScalarCFunc[] sourceFunc, Complex[] param, BVPType tp, Complex[] RHS) {
		assembleGlobalVector(sourceFunc, param, FES1L21::elementCSource, tp, RHS);
	}

	public static void elementSource(int[] ele, double[][] xx, double[][] vf, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xx, 0);
		double L = Math.abs(x[1] - x[0]);
		double[] f = vf[0];
		if (tp == BVPType.COMMON) {
			eleF[0] = (2 * f[0] + f[1]) * L / 6;
			eleF[1] = (f[0] + 2 * f[1]) * L / 6;
		} else if (tp == BVPType.AXISYMMETRIC) {
			eleF[0] = (2 * x[0] * f[0] + MVO.sum(x) * MVO.sum(f)) * L / 12.0;
			eleF[1] = (2 * x[1] * f[1] + MVO.sum(x) * MVO.sum(f)) * L / 12.0;
		} else if (tp == BVPType.SPHERICAL) {
			eleF[0] = (2 * x[0] * x[0] * f[0] + (x[0] * x[0] + x[1] * x[1]) * MVO.sum(f)) * L / 12.0;
			eleF[1] = (2 * x[1] * x[1] * f[1] + (x[0] * x[0] + x[1] * x[1]) * MVO.sum(f)) * L / 12.0;
		}
	}

	public static void elementCSource(int[] ele, double[][] xx, Complex[][] vf, BVPType tp, Complex[] eleF) {
		double[] x = MVO.col(xx, 0);
		double L = Math.abs(x[1] - x[0]);
		Complex[] f = vf[0];
		Complex sumF = Complex.sum(f);
		if (tp == BVPType.COMMON) {
			eleF[0].real = (2 * f[0].real + f[1].real) * L / 6;
			eleF[1].real = (f[0].real + 2 * f[1].real) * L / 6;
			eleF[0].image = (2 * f[0].image + f[1].image) * L / 6;
			eleF[1].image = (f[0].image + 2 * f[1].image) * L / 6;
		}
		else if (tp == BVPType.AXISYMMETRIC) {
			eleF[0].real = (2 * x[0] * f[0].real + MVO.sum(x) * sumF.real) * L / 12.0;
			eleF[1].real = (2 * x[1] * f[1].real + MVO.sum(x) * sumF.real) * L / 12.0;
			eleF[0].image = (2 * x[0] * f[0].image + MVO.sum(x) * sumF.image) * L / 12.0;
			eleF[1].image = (2 * x[1] * f[1].image + MVO.sum(x) * sumF.image) * L / 12.0;
		} else if (tp == BVPType.SPHERICAL) {
			eleF[0].real = (2 * x[0] * x[0] * f[0].real + (x[0] * x[0] + x[1] * x[1]) * sumF.real) * L / 12.0;
			eleF[1].real = (2 * x[1] * x[1] * f[1].real + (x[0] * x[0] + x[1] * x[1]) * sumF.real) * L / 12.0;
			eleF[0].image = (2 * x[0] * x[0] * f[0].image + (x[0] * x[0] + x[1] * x[1]) * sumF.image) * L / 12.0;
			eleF[1].image = (2 * x[1] * x[1] * f[1].image + (x[0] * x[0] + x[1] * x[1]) * sumF.image) * L / 12.0;
		}
	}
	
	public void assembleFlux(ScalarFunc[] coeffFunc, double[] param, BVPType tp, double[] RHS,
							 int... bdLabel) {
		assembleBoundaryVector(coeffFunc, param, FES1L21::elementFlux, tp, RHS, bdLabel);
	}
	
	public void assembleFlux(double[] constCoef, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(constCoef, FES1L21::elementFlux, tp, RHS, bdLabel);
	}

	public void assembleFlux(double[][] vecCoef, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(vecCoef, FES1L21::elementFlux, tp, RHS, bdLabel);
	}
	
	public static void elementFlux(int[] ele, double[][] coord, double[][] g, BVPType tp, double[] k) {
		double x = coord[0][0];
		switch (tp) {
		case COMMON:
			k[0] = g[0][0];
			break;
		case AXIS_PLANESTRAIN:
		case AXIS_PLANESTRESS:
		case AXISYMMETRIC:
			k[0] = g[0][0] * x;
			break;
		case SPHERICAL:
			k[0] = g[0][0] * x * x;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public void assembleCMass(ScalarCFunc[] coefFunc, Complex[] param, BVPType tp, SCMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES1L21::elementCMass, tp, A);
	}
	
	public static void elementMass(int[] ele, double[][] coord, double[][] coef, BVPType tp, double[][] k) {
	    double[] coefM = coef[0];
	    double[] x = MVO.col(coord, 0);
	    double L = Math.abs(x[0] - x[1]);
	    MVO.fill(k, 0.0);
		double xi, wxarea, c;
		double detJ;
	    double[] h = {0.0, 0.0};
	    double[] dXi = {0.0, 0.0};
		double[] coef0, coefX, coefX2, coefX3, coefU, coefU2, coefU3, coefdU, coefdU2, coefdU3, u;
	    switch(tp) {
	    case CONSIST_MASS:
	    	L /= 12.0;
	    	k[0][0] = (3 * coefM[0] + coefM[1]) * L;
	    	k[0][1] = k[1][0] = (coefM[0] + coefM[1]) * L;
	        k[1][1] = (coefM[0] + 3 * coefM[1]) * L;
	        break;
	    case LUMP_MASS:
	    	L /= 6.0;
	    	k[0][0] = (2 * coefM[0] + coefM[1]) * L;
	    	k[1][1] = (coefM[0] + 2 * coefM[1]) * L;
	    	break;
	    case AXIS_CONSIST_MASS:
	    	L /= 12.0;
	    	k[0][0] = (3 * coefM[0] * x[0] + coefM[1] * x[1]) * L;
	    	k[0][1] = k[1][0] = (coefM[0] * x[0] + coefM[1] * x[1]) * L;
	        k[1][1] = (coefM[0] * x[0] + 3 * coefM[1] * x[1]) * L;
	        break;
	    case AXIS_LUMP_MASS:
	        L /= 6.0;
	        k[0][0] = (2 * coefM[0] * x[0] + coefM[1] * x[1]) * L;
	        k[1][1] = (coefM[0] * x[0] + 2 * coefM[1] * x[1]) * L;
	        break;
	    case SPHE_CONSIST_MASS:
	    	L /= 12.0;
	    	k[0][0] = (3 * coefM[0] * x[0] * x[0] + coefM[1] * x[1] * x[1]) * L;
	    	k[0][1] = k[1][0] = (coefM[0] * x[0] * x[0] + coefM[1] * x[1] * x[1]) * L;
	        k[1][1] = (coefM[0] * x[0] * x[0] + 3 * coefM[1] * x[1] * x[1]) * L;
	        break;
	    case SPHE_LUMP_MASS:
	        L /= 6.0;
	        k[0][0] = (2 * coefM[0] * x[0] * x[0] + coefM[1] * x[1] * x[1]) * L;
	        k[1][1] = (coefM[0] * x[0] * x[0] + 2 * coefM[1] * x[1] * x[1]) * L;
	        break;
	    case NONLINEAR_2:
			//系数形式为 a + b x + c x^2 + d x^3 + e u + f u^2 + g u^3 + h du/dx + k (du/dx)^2 + l (du/dx)^3 的非线性切线刚度迭代矩阵
			coefU   = coef[4];
			coefU2  = coef[5];
			coefU3  = coef[6];
			coefdU  = coef[7];
			coefdU2 = coef[8];
			coefdU3 = coef[9];
			u       = coef[10];
//			double dU = (u[1] - u[0]) / L;
//			double avCoef1 = ( 0.5 * (coefdU[0]  + coefdU[1])  * dU
//					+ 1.0 * (coefdU2[0] + coefdU2[1]) * dU * dU
//					+ 1.5 * (coefdU3[0] + coefdU3[1]) * dU * dU * dU ) / L;
//			eleMatrix[0][0] = avCoef1;
//			eleMatrix[0][1] = - avCoef1;
//			eleMatrix[1][0] = - avCoef1;
//			eleMatrix[1][1] = avCoef1;
//
//			double avCoef21 = dU * (2 * coefU[0] + coefU[1] + MVO.sum(coefU2) * MVO.sum(u) + 2 * coefU2[0] * u[0]) / 6;
//			double avCoef22 = dU * (2 * coefU[1] + coefU[0] + MVO.sum(coefU2) * MVO.sum(u) + 2 * coefU2[1] * u[1]) / 6;
//			avCoef21 +=  (12 * coefU3[0] * u[0] * u[0] + 3 * coefU3[1] * u[0] * u[0]
//					+ 6 * coefU3[0] * u[0] * u[1] + 4 * coefU3[1] * u[0] * u[1]
//					+ 2 * coefU3[0] * u[1] * u[1] + 3 * coefU3[1] * u[1] * u[1]) / 20;
//			avCoef22 += (3 * coefU3[0] * u[0] * u[0] + 2 * coefU3[1] * u[0] * u[0]
//					+ 4 * coefU3[0] * u[0] * u[1] + 6 * coefU3[1] * u[0] * u[1]
//					+ 3 * coefU3[0] * u[1] * u[1] + 12 * coefU3[1] * u[1] * u[1]) / 20;
//			eleMatrix[0][0] += -avCoef21 - avCoef22;
//			eleMatrix[0][1] += -avCoef21 - avCoef22;
//			eleMatrix[1][0] += avCoef21 + avCoef22;
//			eleMatrix[1][1] += avCoef21 + avCoef22;
			for(int i = 0; i < 5; i++) {
				xi = Constants.GP5[i];
				detJ = CG.ISOMap1D(CG::L2ShapeFunction, x, xi, h, dXi);
				double c1 = 0;
				double uh = MVO.dot_product(u, h);
				c1 += MVO.dot_product(coefU, h) * uh;
				c1 += 2 * MVO.dot_product(coefU2, h) * uh * uh;
				c1 += 3 * MVO.dot_product(coefU3, h) * uh * uh * uh;

				double c2 = 0.0;
				double dxUh = MVO.dot_product(u, dXi);
				c2 += MVO.dot_product(coefdU, h) * uh;
				c2 += 2 * MVO.dot_product(coefdU2, h) * dxUh * uh;
				c2 += 3 * MVO.dot_product(coefdU3, h) * dxUh * dxUh * uh;
				wxarea = Constants.GW5[i] * detJ;
				for(int si = 0; si < 2; si++) {
					for(int sj = 0; sj < 2; sj++) {
						k[si][sj] += h[si] * h[sj] * wxarea * c1 + h[si] * dXi[sj] * wxarea * c2;
					}
				}
			}
		case NONLINEAR_1:
			//系数形式为 a + b x + c x^2 + d x^3 + e u + f u^2 + g u^3 + h du/dx + k (du/dx)^2 + l (du/dx)^3 的非线性直接迭代矩阵
			coef0   = coef[0];
			coefX   = coef[1];
			coefX2  = coef[2];
			coefX3  = coef[3];
			coefU   = coef[4];
			coefU2  = coef[5];
			coefU3  = coef[6];
			coefdU  = coef[7];
			coefdU2 = coef[8];
			coefdU3 = coef[9];
			u       = coef[10];
//			double avA = MVO.average(coef0) / L;
//			double avBx = (2 * coefX[0] * x[0] + coefX[0] * x[1] + coefX[1] * x[0] + 2 * coefX[1] * x[1]) / (6 * L);
//			double avCx2 = 0.25 * (coefX2[0] * x[0] * x[0] + coefX2[1] * x[1] * x[1]) / L
//					+ (coefX2[0] + coefX2[1]) * x[0] * x[1] / (6 * L)
//					+ (coefX2[0] * x[1] * x[1] + coefX2[1] * x[0] * x[0]) / (12 * L);
//			double avDx3 = ( (4 * coefX3[0] + coefX3[1]) * x[0] * x[0] * x[0]
//					+ (3 * coefX3[0] + 2 * coefX3[1]) * x[0] * x[0] * x[1]
//					+ (2 * coefX3[0] + 3 * coefX3[1]) * x[0] * x[1] * x[1]
//					+ (coefX3[0] + 4 * coefX3[1]) * x[1] * x[1] * x[1]) / (20 * L);
//			double avEu = (2 * coefU[0] * u[0] + coefU[0] * u[1] + coefU[1] * u[0] + 2 * coefU[1] * u[1]) / (6 * L);
//			double avFu2 = 0.25 * (coefU2[0] * u[0] * u[0] + coefU2[1] * u[1] * u[1]) / L
//					+ (coefU2[0] + coefU2[1]) * u[0] * u[1] / (6 * L)
//					+ (coefU2[0] * u[1] * u[1] + coefU2[1] * u[0] * u[0]) / (12 * L);
//			double avGu3 = ((4 * coefU3[0] + coefU3[1]) * u[0] * u[0] * u[0]
//					+ (3 * coefU3[0] + 2 * coefU3[1]) * u[0] * u[0] * u[1]
//					+ (2 * coefU3[0] + 3 * coefU3[1]) * u[0] * u[1] * u[1]
//					+ (coefU3[0] + 4 * coefU3[1]) * u[1] * u[1] * u[1]) / (20 * L);
//			dxU = (x[1] - x[0]) / L;
//			double avHDu  = 0.5 * (coefdU[0] + coefdU[1]) * dxU / L;
//			double avKDu2 = 0.5 * (coefdU2[0] + coefdU2[1]) * dxU * dxU / L;
//			double avLDu3 = 0.5 * (coefdU3[0] + coefdU3[1]) * dxU * dxU * dxU / L;
//			double totalCoef = avA + avBx + avCx2 + avDx3 + avEu + avFu2 + avGu3 + avHDu + avKDu2 + avLDu3;
//			eleMatrix[0][0] += totalCoef;
//			eleMatrix[0][1] += - totalCoef;
//			eleMatrix[1][0] += - totalCoef;
//			eleMatrix[1][1] += totalCoef;
			for(int i = 0; i < 5; i++) {
				xi = Constants.GP5[i];
				detJ = CG.ISOMap1D(CG::L2ShapeFunction, x, xi, h, dXi);
				c = MVO.dot_product(coef0, h);
				double xh = MVO.dot_product(x, h);
				c += MVO.dot_product(coefX, h) * xh;
				c += MVO.dot_product(coefX2, h) * xh * xh;
				c += MVO.dot_product(coefX3, h) * xh * xh * xh;
				double uh = MVO.dot_product(u, h);
				c += MVO.dot_product(coefU, h) * uh;
				c += MVO.dot_product(coefU2, h) * uh * uh;
				c += MVO.dot_product(coefU3, h) * uh * uh * uh;
				double dxUh = MVO.dot_product(u, dXi);
				c += MVO.dot_product(coefdU, h) * dxUh;
				c += MVO.dot_product(coefdU2, h) * dxUh * dxUh;
				c += MVO.dot_product(coefdU3, h) * dxUh * dxUh * dxUh;
				wxarea = Constants.GW5[i] * detJ * c;
				for(int si = 0; si < 2; si++) {
					for(int sj = 0; sj < 2; sj++) {
						k[si][sj] += h[si] * h[sj] * wxarea;
					}
				}
			}
			break;
	    default:
	    	throw new IllegalArgumentException("Wrong BVPType!");
	    }
	    //System.out.println(MVO.toString(k));
	    //System.exit(0);
	}

	public static void elementCMass(int[] ele, double[][] coord, Complex[][] g, BVPType tp, Complex[][] k) {
		Complex[] coefM = new Complex[2];
		for(int i = 0; i < 2; i++) {
			coefM[i] = new Complex(g[0][i]);
		}
		double[] x = MVO.col(coord, 0);
		double L = Math.abs(x[0] - x[1]);
		switch (tp) {
			case CONSIST_MASS -> {
				L /= 12.0;
				k[0][0].real = (3 * coefM[0].real + coefM[1].real) * L;
				k[0][1].real = k[1][0].real = (coefM[0].real + coefM[1].real) * L;
				k[1][1].real = (coefM[0].real + 3 * coefM[1].real) * L;
				k[0][0].image = (3 * coefM[0].image + coefM[1].image) * L;
				k[0][1].image = k[1][0].image = (coefM[0].image + coefM[1].image) * L;
				k[1][1].image = (coefM[0].image + 3 * coefM[1].image) * L;
			}
			case LUMP_MASS -> {
				L /= 6.0;
				k[0][0].real = (2 * coefM[0].real + coefM[1].real) * L;
				k[1][1].real = (coefM[0].real + 2 * coefM[1].real) * L;
				k[0][1].real = k[1][0].real = k[0][1].image = k[1][0].image = 0.0;
				k[0][0].image = (2 * coefM[0].image + coefM[1].image) * L;
				k[1][1].image = (coefM[0].image + 2 * coefM[1].image) * L;
			}
			case AXIS_CONSIST_MASS -> {
				L /= 12.0;
				k[0][0].real = (3 * coefM[0].real * x[0] + coefM[1].real * x[1]) * L;
				k[0][1].real = k[1][0].real = (coefM[0].real * x[0] + coefM[1].real * x[1]) * L;
				k[1][1].real = (coefM[0].real * x[0] + 3 * coefM[1].real * x[1]) * L;
				k[0][0].image = (3 * coefM[0].image * x[0] + coefM[1].image * x[1]) * L;
				k[0][1].image = k[1][0].image = (coefM[0].image * x[0] + coefM[1].image * x[1]) * L;
				k[1][1].image = (coefM[0].image * x[0] + 3 * coefM[1].image * x[1]) * L;
			}
			case AXIS_LUMP_MASS -> {
				L /= 6.0;
				k[0][0].real = (2 * coefM[0].real * x[0] + coefM[1].real * x[1]) * L;
				k[1][1].real = (coefM[0].real * x[0] + 2 * coefM[1].real * x[1]) * L;
				k[0][0].image = (2 * coefM[0].image * x[0] + coefM[1].image * x[1]) * L;
				k[1][1].image = (coefM[0].image * x[0] + 2 * coefM[1].image * x[1]) * L;
				k[0][1].real = k[1][0].real = k[0][1].image = k[1][0].image = 0.0;
			}
			case SPHE_CONSIST_MASS -> {
				L /= 12.0;
				k[0][0].real = (3 * coefM[0].real * x[0] * x[0] + coefM[1].real * x[1] * x[1]) * L;
				k[0][1].real = k[1][0].real = (coefM[0].real * x[0] * x[0] + coefM[1].real * x[1] * x[1]) * L;
				k[1][1].real = (coefM[0].real * x[0] * x[0] + 3 * coefM[1].real * x[1] * x[1]) * L;
				k[0][0].image = (3 * coefM[0].image * x[0] * x[0] + coefM[1].image * x[1] * x[1]) * L;
				k[0][1].image = k[1][0].image = (coefM[0].image * x[0] * x[0] + coefM[1].image * x[1] * x[1]) * L;
				k[1][1].image = (coefM[0].image * x[0] * x[0] + 3 * coefM[1].image * x[1] * x[1]) * L;
			}
			case SPHE_LUMP_MASS -> {
				L /= 6.0;
				k[0][0].real = (2 * coefM[0].real * x[0] * x[0] + coefM[1].real * x[1] * x[1]) * L;
				k[1][1].real = (coefM[0].real * x[0] * x[0] + 2 * coefM[1].real * x[1] * x[1]) * L;
				k[0][0].image = (2 * coefM[0].image * x[0] * x[0] + coefM[1].image * x[1] * x[1]) * L;
				k[1][1].image = (coefM[0].image * x[0] * x[0] + 2 * coefM[1].image * x[1] * x[1]) * L;
				k[0][1].real = k[1][0].real = k[0][1].image = k[1][0].image = 0.0;
			}
			default -> throw new IllegalArgumentException("Wrong BVPType!");
		}
		//System.out.println(MVO.toString(k));
		//System.exit(0);
	}

	public void assembleConvection(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES1L21::elementConvection, tp, A);
	}

	public void assembleConvection(double[] constCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(constCoef, FES1L21::elementConvection, tp, A);
	}

	public void assembleConvection(double[][] vecCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(vecCoef, FES1L21::elementConvection, tp, A);
	}

	public static void elementConvection(int[] ele, double[][] xx, double[][] kxx, BVPType tp, double[][] eleMatrix) {
//		double[] x = MVO.col(xx, 0);
//		double L = Math.abs(x[0] - x[1]);
		double coef1 = 0.0;
		double coef2 = 0.0;
		double sumB = MVO.sum(kxx[0]);
		if (tp == BVPType.COMMON) {
			coef1 = (sumB + kxx[0][0]) / 6;
			coef2 = (sumB + kxx[0][1]) / 6;
		} else if (tp == BVPType.AXISYMMETRIC){
			//coef = (MVO.sum(kxx[0]) * MVO.sum(x) + MVO.dot_product(kxx[0], x)) / (6.0 * L);
		} else {
			throw new IllegalArgumentException("Wrong BVPType!");
		}
		eleMatrix[0][0] = -coef1;
		eleMatrix[0][1] = coef1;
		eleMatrix[1][0] = -coef2;
		eleMatrix[1][1] = coef2;
	}
	
	public void assembleBoundaryMass(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A, int... bdLabel) {
		assembleBoundaryMatrix(coefFunc, param, FES1L21::elementConvectSurface, tp, A, bdLabel);
	}
	
	public void assembleBoundaryMass(double[] constConvect, BVPType tp, SRMatrix A, int... bdLabel) {
		assembleBoundaryMatrix(constConvect, FES1L21::elementConvectSurface, tp, A, bdLabel);
	}
	
	public void assembleBoundaryMass(double[][] vectorConvect, BVPType tp, SRMatrix A, int... bdLabel) {
		assembleBoundaryMatrix(vectorConvect, FES1L21::elementConvectSurface, tp, A, bdLabel);
	}
	
	public static void elementConvectSurface(int[] ele, double[][] coord, double[][] h, BVPType tp, double[][] eleMatrix) {
		double x = coord[0][0];
		switch (tp) {
		case CONSIST_MASS:
			eleMatrix[0][0] = h[0][0];
			break;
		case AXIS_CONSIST_MASS:
		case AXIS_LUMP_MASS:
			eleMatrix[0][0] = h[0][0] * x;
			break;
		case SPHE_CONSIST_MASS:
		case SPHE_LUMP_MASS:
			eleMatrix[0][0] = h[0][0] * x * x;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
	public void assembleHeatLoad(ScalarFunc[] coef, double[] param, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, param, vecCoef, FES1L21::elementHeatLoad, tp, RHS);
	}

	public void assembleHeatLoad(double[] coef, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, vecCoef, FES1L21::elementHeatLoad, tp, RHS);
	}

	public static void elementHeatLoad(int[] ele, double[][] coord, double[][] coef, BVPType tp, double[] eleLoad) {
	    double[] beta = coef[0];
	    double[] heat = coef[1];
	    double bh = 0.5 * (beta[0] * heat[0] + beta[1] * heat[1]);
	    double[] x = MVO.col(coord, 0);
	    double L = Math.abs(x[1] - x[0]);
	    double L0 = 0.50 * MVO.sum(x);
	    double[] dx = {-1.0/L, 1.0/L};
	    if (tp == BVPType.COMMON) {
	        eleLoad[0] = -bh;
	        eleLoad[1] = bh;
	    } else if(tp == BVPType.AXISYMMETRIC) {
	        eleLoad[0] = bh * (dx[0] + 0.50 / L0) * L * L0;
	        eleLoad[1] = bh * (dx[1] + 0.50 / L0) * L * L0;
	    } else if (tp == BVPType.SPHERICAL) { //这个还有问题
	        eleLoad[0] = bh * (dx[0] + 0.50 / L0) * L * L0;
	        eleLoad[1] = bh * (dx[1] + 0.50 / L0) * L * L0;	    	
	    } else {
	    	throw new IllegalArgumentException("Wrong BVPType!");
	    }
	}
	
	
	public void computeGradient(double[] u, double[][] gradient) {

	    double[] count = new double[ndof];
	    double[][] coord = new double[2][1];
	    double[] eleU = new double[2];
	    MVO.fill(gradient, 0.0);
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	double[] x = MVO.col(coord, 0);
	    	eleU[0] = u[ele[0]];
	    	eleU[1] = u[ele[1]];
	    	gradient[0][ele[0]] += eleU[1] - eleU[0];
	    	gradient[0][ele[1]] += eleU[1] - eleU[0];
	    	count[ele[0]] += x[1] - x[0];
	    	count[ele[1]] += x[1] - x[0];
	    }
	    for(int i = 0; i < ndof; i++) {
	    	gradient[0][i] /= count[i];
	    }
	    
	    double[] xx = { mesh.getNode(0)[0], mesh.getNode(1)[0], mesh.getNode(2)[0]};
	    //�ֱ����߽��ϵĵ���dU(1)��dU(N)
	    double dx1 = xx[1] - xx[0]; 
	    double dx2 = xx[2] - xx[1];
	    gradient[0][0] = ((u[1] - u[0]) * (dx1 + dx2) * (dx1 + dx2) 
	    		        - (u[2] - u[0]) * dx1 * dx1) / (dx1 * dx2 * (dx1 + dx2));
	    
	    int N = ndof;
	    xx[0] = mesh.getNode(N - 3)[0];
	    xx[1] = mesh.getNode(N - 2)[0];
	    xx[2] = mesh.getNode(N - 1)[0];
	    //dx1 = x(N) - x(N-1); dx2 = x(N - 1) - x(N - 2)
	    dx1 = xx[2] - xx[1]; dx2 = xx[1] - xx[0];
	    gradient[0][N - 1] = -((u[N - 2] - u[N - 1]) * (dx1 + dx2)*(dx1 + dx2)  
	    		             - (u[N - 3] - u[N - 1]) * dx1*dx1) / (dx1 * dx2 * (dx1 + dx2));
	}

	public Pair<double[], double[]> computeElementGradient(double[] u) {
		double[] count = new double[ndof];
		double[][] coord = new double[2][1];
		double[] eleU = new double[2];
		int n = mesh.getNt();
		double[] xCenter = new double[n];
		double[] gradient = new double[n];
		for(int i = 0; i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			double[] x = MVO.col(coord, 0);
			eleU[0] = u[ele[0]];
			eleU[1] = u[ele[1]];
			gradient[i] = (eleU[1] - eleU[0]) / (x[1] - x[0]);
			xCenter[i] = 0.5 * (x[1] + x[0]);
		}
		return new Pair<>(xCenter, gradient);
	}

	public double[] computeGradient(double[] u) {
		double[] coordX = mesh.getCoords()[0];
		Pair<double[], double[]> eleInfo = computeElementGradient(u);
		return PDApproximation.Interpolation1D(eleInfo.getFirst(), eleInfo.getSecond(), coordX, 2);
	}

	public double computeL2Norm(double[] sol) {
	    double sum = 0.0;
	    double[][] coord = new double[2][1];
	    double[] u = new double[2];
	    for(int i = 0, nt = mesh.getNt(); i < nt; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	double length = Math.abs(coord[0][0] - coord[1][0]);
	    	u[0] = sol[ele[0]];
	    	u[1] = sol[ele[1]];
	    	sum += length * (u[0] * u[0] + u[1] * u[1] + u[0] * u[1]) / 3.0;
	    }
	    return Math.sqrt(sum);
	}

	public double computeH1SemiNorm(double[] sol) {
	    double sum = 0.0;
	    double[][] coord = new double[2][1];
	    double[] u = new double[2];
	    double length, du;
	    for(int i = 0, nt = mesh.getNt(); i < nt; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	double[] x = MVO.col(coord, 0);
	    	length = Math.abs(x[1] - x[0]);
	    	u[0] = sol[ele[0]];
	    	u[1] = sol[ele[1]];
	        du = (u[1] - u[0])/(x[1] - x[0]);
	        sum += length * du * du;
	    }
	    return Math.sqrt(sum);
	}
	
	public double integrate(double[] sol) {
		double result = 0.0;
		double[][] coord = new double[2][1];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
	    	double[] x = MVO.col(coord, 0);
	    	double length = Math.abs(x[1] - x[0]);
			result += 0.5 * (sol[ele[0]] + sol[ele[1]]) * length;
		}
		return result;
	}

	public void applyABC(SCMatrix A, double k0, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		for(int idx: bcSet) {
			A.addElement(idx, idx, new Complex(0, k0));
		}
	}

	public double[] computeReactionForce(double[][] vecCoef, BVPType tp, double[] solution) {
		return computeReaction(vecCoef, tp, FES1L21::elementStiff, solution);
	}

	public double[] computeReactionForce(ScalarFunc[] coeffunc, double[] param, BVPType tp, double[] solution) {
		return computeReaction(coeffunc, param, tp, FES1L21::elementStiff, solution);
	}

	public double[] computeReactionForce(Function<Double, Double> forceFunc, BVPType tp, double[] solution) {
		double[] solDof = new double[dofPerElement];
		double[] result = new double[solution.length];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			int[] dof = getDof(ele);
			mesh.getCoordPerElement(ele, coordRegion);
			MVO.getSlice(solution, dof, solDof);
			double strain = (solDof[0] - solDof[1]) / (coordRegion[0][0] - coordRegion[1][0]);
			double force = forceFunc.apply(strain);
			double[] reaction = {-force, force};
			assembleVector(reaction, dof, result);
		}
		return result;
	}

	public Pair<double[], double[]> computeElementStress(ScalarFunc elasticCoef, double[] param, double[] u) {
		double[] count = new double[ndof];
		double[][] coord = new double[2][1];
		double[] eleU = new double[2];
		int n = mesh.getNt();
		double[] xCenter = new double[n];
		double[] stress = new double[n];
		for(int i = 0; i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			double[] x = MVO.col(coord, 0);
			eleU[0] = u[ele[0]];
			eleU[1] = u[ele[1]];
			double E = elasticCoef.action(x, mesh.getElementLabel(i), param);
			stress[i] = E * (eleU[1] - eleU[0]) / (x[1] - x[0]);
			xCenter[i] = 0.5 * (x[1] + x[0]);
		}
		return new Pair<>(xCenter, stress);
	}

	public double[] computeElastoplasticActionForce(ScalarFunc[] elastoPlastoCoef, double[] param, BVPType tp, double[] u) {
		double[][] coef = new double[elastoPlastoCoef.length][nPerEle];
		double[] result = new double[u.length];
		double[][] coord = new double[2][1];
		double[] eleU = new double[2];
		int n = mesh.getNt();
		double force = 0.0;
		for(int i = 0; i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			getCoefFromFunc(elastoPlastoCoef, coordRegion, mesh.getElementLabel(i), param, coef);
			int[] dof = getDof(ele);
			eleU[0] = u[ele[0]];
			eleU[1] = u[ele[1]];
			double E = MVO.average(coef[0]);
			double A = MVO.average(coef[1]);
			double yield = MVO.average(coef[2]);
			double Hards = MVO.average(coef[3]);
			double[] x = MVO.col(coord, 0);
			double dStrain = (eleU[1] - eleU[0]) / (x[1] - x[0]);
			double stlin = E * dStrain; //线弹性应力增量
			double stcur = elementStress[i] + stlin; //在弹性加载下当前步的应力
			double preys = yield + Hards * Math.abs(elementPlasticStrain[i]); //之前的屈服应力
			if(Math.abs(elementStress[i]) < preys) { //如果当前单元在之前步还没有发生屈服， 则判断当前步单元是否发生了屈服
				if(Math.abs(stcur) < preys) { //如果当前应力仍小于屈服应力, 则不需要作处理, 将当前应力赋予单元应力
					elementStress[i] = stcur;
				} else {  //发生了屈服
					double R = (Math.abs(stcur) - preys) / Math.abs(stlin);
					if (R < 0.0 || elementStress[i] * stlin < 0.0) { //是否是卸载过程
						elementStress[i] = stcur;
					} else{
						//初始屈服时更新应力与塑性应变
						elementStress[i] += (1 - R) * stlin + R * E * (1.0 - E / (E + Hards)) * dStrain;
						elementPlasticStrain[i] += R * dStrain * E / (E + Hards);
					}
				}
			} else { //对于已经发生屈服的单元
				if( elementStress[i] * stlin < 0.0) { //判断是否是卸载过程
					elementStress[i] = stcur;
				} else {
					double R = 1.0;
					elementStress[i] += (1 - R) * stlin + R * E * (1.0 - E / (E + Hards)) * dStrain;
					elementPlasticStrain[i] += R * dStrain * E / (E + Hards);
				}
			}
			//当前单元力为应力乘以面积
			force = elementStress[i] * A;
			double sign = Math.signum(x[1] - x[0]);
			double[] reaction = {-sign * force, sign * force};
			assembleVector(reaction, dof, result);
		}
		return result;
	}

	double[] elementElasticStrain;
	double[] elementStress;
	double[] elementPlasticStrain;

}
